function VignPlus(X,Y)
    error if #X ne #Y, "Arguments must have the same length.";
    return &*[ C2S(((S2C(X[i])+S2C(Y[i])-130) mod 26)+65) : i in [1..#X] ];
end function;

function VignDiff(X,Y)
    error if #X ne #Y, "Arguments must have the same length.";
    return &*[ C2S(((S2C(X[i])-S2C(Y[i])) mod 26)+65) : i in [1..#X] ];
end function;

function DigraphDifferentialProbabilities1(DD,PT)
    assert #DD eq 2;
    AZ := {@ CodeToString(64+i) : i in [1..26] @};
    r := Index(AZ,DD[1]); s := Index(AZ,DD[2]);
    FD2 := DigraphFrequencyDistribution(PT);
    DD2 := Parent(FD2)!0;
    FDD := RealField()!0;
    for i1, j1 in [1..26] do
	i2 := ((i1-r) mod 26) + 1;
	j2 := ((j1-s) mod 26) + 1;
	F2 := FD2[i1,j1] * FD2[i2,j2];
	DD2[i1,j1] +:= F2;
	FDD +:= F2;
    end for;
    return (1/FDD)*DD2;
end function;

function DigraphDifferentialProbabilities2(D1,D2,PT)
    assert #D1 eq 2 and #D2 eq 2;
    AZ := {@ CodeToString(64+i) : i in [1..26] @};
    r1 := Index(AZ,D1[1]); s1 := Index(AZ,D1[2]);
    r2 := Index(AZ,D2[1]); s2 := Index(AZ,D2[2]);
    FDD := RealField()!0;
    DD2 := MatrixAlgebra(RealField(),26)!0;
    F2D := DigraphFrequencyDistribution(PT);
    for i1, j1 in [1..26] do
	i2 := ((i1-r1) mod 26) + 1;
	j2 := ((j1-s1) mod 26) + 1;
	i3 := ((i1-r2) mod 26) + 1;
	j3 := ((j1-s2) mod 26) + 1;
	F3 := F2D[i1,j1] * F2D[i2,j2] * F2D[i3,j3];
	DD2[i1,j1] +:= F3;
	FDD +:= F3;
    end for;
    error if FDD eq 0,
	"Arguments 1 and 2 must have nonzero probabilities.";
    return (1/FDD)*DD2;
end function;

function TrigraphDifferentialProbabilities1(DD,PT)
    assert #DD eq 3;
    AZ := {@ CodeToString(64+i) : i in [1..26] @};
    r := Index(AZ,DD[1]);
    s := Index(AZ,DD[2]);
    t := Index(AZ,DD[2]);
    AAA, FD3 := PolygraphFrequencyDistribution(PT,3,RealField()!0);
    DD3 := [ RealField()!0 : i in [1..26^3] ];
    FDD := RealField()!0;
    for i, j, k in [1..26] do
	l := ((i-r) mod 26) + 1;
	m := ((j-s) mod 26) + 1;
	n := ((k-t) mod 26) + 1;
	XYZ := AZ[i]*AZ[j]*AZ[k]; 
	WQR := AZ[l]*AZ[m]*AZ[n]; 
	x := Index(AAA,XYZ);
	y := Index(AAA,WQR);
	if x ne 0 and y ne 0 then
	    F3 := FD3[x]*FD3[y];
	    DD3[26^2*(i-1)+26*(j-1)+k] +:= F3;
	    FDD +:= F3;
	end if;
    end for;
    error if FDD eq 0, "Argument 1 must have nonzero probability.";
    return [ (1/FDD)*f : f in DD3 ];
end function;

function TrigraphDifferentialProbabilities2(D1,D2,PT)
    assert #D1 eq 3 and #D2 eq 3;
    AZ := {@ CodeToString(64+i) : i in [1..26] @};
    r1 := Index(AZ,D1[1]); s1 := Index(AZ,D1[2]); t1 := Index(AZ,D1[3]);
    r2 := Index(AZ,D2[1]); s2 := Index(AZ,D2[2]); t2 := Index(AZ,D2[3]);
    FDD := RealField()!0;
    DD3 := [ RealField()!0 : i in [1..26^3] ];
    AAA, FD3 := PolygraphFrequencyDistribution(PT,3,RealField()!0);
    for i1, j1, k1 in [1..26] do
	///////////////////////////
	i2 := ((i1-r1) mod 26) + 1;
	j2 := ((j1-s1) mod 26) + 1;
	k2 := ((k1-t1) mod 26) + 1;
	///////////////////////////
	i3 := ((i1-r2) mod 26) + 1;
	j3 := ((j1-s2) mod 26) + 1;
	k3 := ((k1-t2) mod 26) + 1;
	///////////////////////////
	XYZ := AZ[i1]*AZ[j1]*AZ[k1]; 
	QRS := AZ[i2]*AZ[j2]*AZ[k2]; 
	TUV := AZ[i3]*AZ[j3]*AZ[k3]; 
	x := Index(AAA,XYZ);
	y := Index(AAA,QRS);
	z := Index(AAA,TUV);
	if x ne 0 and y ne 0 and z ne 0 then
	    F3 := FD3[x] * FD3[y] * FD3[z];
	    FDD +:= F3;
	    DD3[26^2*(i1-1)+26*(j1-1)+k1] +:= F3;
	end if;
    end for;
    error if FDD eq 0, "Argument 1 must have nonzero probability.";
    return [ (1/FDD)*f : f in DD3 ];
end function;

////////////////////////////////////////////////////////////////////////////////

C2S := CodeToString;
S2C := StringToCode;

AZ := {@ CodeToString(64+i) : i in [1..26] @};

PT := StripEncoding(Read("txt/PoeEdgarAllan/blackcat.txt"));
PT1 := "TOBEORNOTTOBE";
PT2 := "THECATINTHEHA";
PT3 := "WHENINTHECOUR";
DD1 := "AHXCOYFBAMKUE"; // PT1-PT2;
DD2 := "XHXRGEUHPRAHN"; // PT1-PT3;
DD3 := "XAAPSGPGPFQNJ"; // PT2-PT3;

/*
n := #PT-1;
AA := {@ X*Y : X, Y in AZ @}
    where AZ := [ C2S(64+i) : i in [1..26] ];
n := 2001;
eps := 1.0/n^2;
D2 := [ RealField() | 0 : i in [1..26^2] ];
tyme := Cputime();
for i in [1..n] do
    for j in [1..n] do
	if i ne j then
	    DD := VignDiff(PT[[i,i+1]],PT[[j,j+1]]);
	    D2[Index(AA,DD)] +:= eps;
	end if;
    end for;
    if i mod 100 eq 0 then
	printf "%5o: %o\n", i, Cputime(tyme);
    end if;
end for;
*/

for i in [1..12] do
    D1 := DD1[[i,i+1]];
    FD2 := DigraphDifferentialProbabilities1(D1,PT);
    printf "%2o:\n", i;
    printf "    %o\n", [ PT1[[i..i+1]], PT2[[i..i+1]] ];
    [ [XY,ZW]
	where ZW := VignDiff(XY,D1)
	where XY := AZ[i]*AZ[j] : i, j in [1..26] | FD2[i,j] ge 0.05 ];
end for;
    
for i in [1..12] do
    D1 := DD1[[i,i+1]];
    D2 := DD2[[i,i+1]];
    FD2 := DigraphDifferentialProbabilities2(D1,D2,PT);
    printf "%2o:\n", i;
    printf "    %o\n", [ PT1[[i..i+1]], PT2[[i..i+1]], PT3[[i..i+1]] ];
    [ [XY,ZW,QR]
	where QR := VignDiff(XY,D2)
	where ZW := VignDiff(XY,D1)
	where XY := AZ[i]*AZ[j] : i, j in [1..26] | FD2[i,j] ge 0.1 ];
end for;
    
for i in [1..11] do
    D1 := DD1[[i..i+2]];
    D2 := DD2[[i..i+2]];
    FD3 := TrigraphDifferentialProbabilities2(D1,D2,PT);
    i0 := Index(AZ,PT1[i]);
    j0 := Index(AZ,PT1[i+1]);
    k0 := Index(AZ,PT1[i+2]);
    printf "%2o:\n", i;
    printf "    %o (%o)\n",
	[ PT1[[i..i+2]], PT2[[i..i+2]], PT3[[i..i+2]] ],
	FD3[26^2*(i0-1)+26*(j0-1)+k0];
    [ [XYZ,QRS,TUV]
	where TUV := VignDiff(XYZ,D2)
	where QRS := VignDiff(XYZ,D1)
	where XYZ := AZ[i]*AZ[j]*AZ[k] :
	i, j, k in [1..26] | FD3[26^2*(i-1)+26*(j-1)+k] ge 0.2 ];
end for;
    
